function [ u ] = evaluateInteriorPoint( x, y, z, mesh, theta )
%EVALUATEINTERIORPOINT Summary of this function goes here
%   Detailed explanation goes here

n = size(x, 1);
nTri = size(mesh.f, 1);

u = zeros(n, 3);

for i = 1 : n
    
    p = [x(i), y(i), z(i)];
    
    for j = 1 : nTri              
        
        q0 = mesh.v(mesh.f(j, 1), :);
        q1 = mesh.v(mesh.f(j, 2), :);
        q2 = mesh.v(mesh.f(j, 3), :);
        
        T = triarea(q0, q1, q2);
              
        qn = mesh.fn(j, :);
        
        v0 = evaluateQuadraturePoint(p, q0);
        v1 = evaluateQuadraturePoint(p, q1);
        v2 = evaluateQuadraturePoint(p, q2);
        
        u(i, :) = u(i, :) + ((v0 + v1 + v2) * T / 3.0 * theta(3*j-2:3*j))';
        
%         u(i, :) = u(i, :) + ((v0 + v1) * domain2D.l(j) / 2.0 * theta(2*j-1:2*j))'; 
    end
    
end

end

